#set working directory, clear workspace
setwd("")
rm(list=ls())

#installing/loading necessary pacakges
#install.packages("R2WinBUGS")
#install.packages("foreign")
#install.packages("mcmcplots")
#install.packages("plyr")
library(coda)
library(foreign)
library(R2WinBUGS)
library(mcmcplots)
library(plyr)

##################################################################################################################################  

#ENGLAND

#reading in matrix of votes and economic variables
dat <- read.csv("england-vote-data-final.csv")
dat <- na.omit(dat)

#defining variable names to match BUGS file
  y <- dat$currentrate + dat$vote
  nvote <- length(y)
  bank <- dat$cbid
  nbank <- nrow(count(dat,vars="cbid"))
  sig <- dat$width.vote
  lag <- dat$currentrate
  ygap <- dat$ygap
  igap <- dat$infgap
  xr <- dat$xr

#reading in the model file
  model.file <- "rcm-corr-uk.txt"

#setting intial values
  inits <- function() { list(a=runif(nbank,-1,1),
                           bsig=runif(nbank,-1,1),
                           blag=runif(nbank,-1,1),
                           bygap=rnorm(1),
                           bigap=rnorm(1),
                           bxr=rnorm(1),
                           
                           sigma.a=runif(1),
                           sigma.bsig=runif(1),
                           rho=runif(1),
                           sigma.blag=runif(1),
                           sigma.bygap=runif(1),
                           sigma.bigap=runif(1),
                           sigma.bxr=runif(1)
                            )}

#running model in WinBUGS
  bugs.model <- bugs(data=c("y","nvote","bank","nbank", "sig","lag","ygap","igap","xr"),
                   inits=inits,
                   parameters.to.save=c("a","bsig","blag","bygap","bigap","bxr",
                                        "mu.a","mu.bsig","mu.blag","mu.bygap","mu.bigap","mu.bxr",
                                        "sigma.a","sigma.bsig","sigma.blag","sigma.bygap","sigma.bigap","sigma.bxr","rho"),
                   model.file= model.file,
                   n.chains=3,
                   n.iter=100000,
                   n.burnin=50000,
                   n.thin=5,
                   debug=TRUE)

#Extracting Coefficients from bugs.model for Table 1:
  bugs.model #Deviance, DIC
  bugs.model$mean$mu.a
  bugs.model$sd$mu.a
  bugs.model$mean$mu.bsig
  bugs.model$sd$mu.bsig
  bugs.model$mean$mu.blag
  bugs.model$sd$mu.blag
  bugs.model$mean$bygap
  bugs.model$sd$bygap
  bugs.model$mean$bigap
  bugs.model$sd$bigap

#Saving Individual Preference Estimates
  b.int <- bugs.model$mean$a

  names <- c("W. Buiter", "E. George", "M. King","A. Budd","C. Goodhart", "J. Vickers", "D. Julius", "I. Plenderleith",
           "D. Clementi", "S. Wadhwani","S.Nickell","C. Allsopp", "C. Bean","P. Tucker","M. Bell", "A. Large", "R. Lambert",
           "R. Lomax", "D. Walton", "J. Gieve", "D. Blanchflower", "T. Besley", "A. Sentence", "S. Dale","P. Fisher",
           "D. Miles", "A. Posen", "M. Weale", "B. Broadbent","I. McCafferty", "M. Carney", "J. Cunliffe", "A. Haldane", 
           "K. Forbes", "N. Shafik", "K. Barker", "M. Saunders", "G. Vlieghe")

  pref.parameters <- cbind(names, b.int)

#Adding Information about Appointing Party and the Level of Uncertainty At the Time of Appointment
  appoint.dat <- read.csv("appointdata.csv")
    
    party.appoint <- appoint.dat$partycode
    appoint.year <- appoint.dat$appoint.year
    appoint.width <- appoint.dat$appoint.width
    
pref.parameters <- cbind(pref.parameters,party.appoint,appoint.year,appoint.width)
pref.parameters <- as.data.frame(pref.parameters)

write.csv(pref.parameters,file="England_Pref.csv")

###################################################################################################################################
  
dat.pref <- read.csv("England_Pref.csv")

#Estimate Appointment Models
  model1 <- lm(dat.pref$b.int~dat.pref$appoint.width)
    summary(model1)

  model2 <- lm(dat.pref$b.int~dat.pref$appoint.width+factor(dat.pref$party.appoint)-1)
    summary(model2)

###################################################################################################################################

    